Quasi-harmonic vs. 



"exact" surface free energies of Al: a systematic study employing 
a new interatomic potential 
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We discuss a computationally efficient classical many-body potential designed to model the Al- 
Al interaction in a wide range of bonding geometries. We show that the potential yields results 
in properties in excellent agreement with experiment and ab initio results for a number of bulk 
and surface properties, among others for surface and step formation energies, and self-diffusion 
barriers. As an application, free energy calculations are performed for the Al (100) surface by 
Monte Carlo thermodynamic integration and the quasi-harmonic approximation. Comparison of 
the latter approximation with the reference Monte Carlo results provides informations on its range 
of applicability to surface problems at high temperatures. 
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I. INTRODUCTION 



Atomistic simulations are playing an increasingly 
prominent role in materials science. From studies of 
crystallisation of clustersa to large-scale simulations of 
fracturea and grain boundary diffusiorfcl, atomistic sim- 
ulations offer a microscopic physical view that cannot 
be obtained from experiment. Predictions resulting from 
this atomic level understanding are proving increasingly 
accurate and useful a 

The effective interatomic interaction potential is the 
key ingredient in all atomistic simulation. The accuracy 
of the potential affects drastically the quality of the sim- 
ulation result, and its functional complexity determines 
the amount of computer time required. U Much research 
effort has therefore beeai devoted to the design of po- 
tential energy functions!! This is especially important in 
classical dynamics which, although quantum mechanical 
simulations have been progressing at a rapid pace in re- 
cent years, remains the most (sometimes, the only) af- 
fordable way to perform very large scale simulations in 
materials science. In this paper, we present a new care- 
fully designed Al-Al interaction model, test its perfor- 
mance, and we apply it to the study of free energies in 
atomic scale simulations. 

The ability to compute free energies is essential to un- 
derstand or predict many physical phenomena, from the 
stability of crystal structures, to the propensity to form 
defects or disorder, and to morphology changes and phase 
transitions. However, the determination of free energies 
from atomic scale computer simulations is a daunting 
task. Approximate methods, mostly based on the har- 
monic vibrational properties of the system, are commonly 
in use to this end. Here, we compare several possible 
versions of the so called quasi-harmonic approximation, 
using as reference accurate simulation using canonical or 
constant-pressure Monte Carlo methods and thermody- 
namic integration, focusing on the specific case of surface 



free energies. The goal is to provide a measure of the 
range of applicability of approximate methods for com- 
plex systems using a reliable Al interaction model. 



II. AL INTERACTION POTENTIAL 

Previously developed interatomic potentialsi~0 for the 
Al-Al interaction have mostly focused on bulk and molec- 
ular properties. In this work, we analyze and generalize 
one of those models with special regard to surface prop- 
erties, aiming as usual at describing Al in as wide a range 
of chemical environments as possible, i.e. ranging from 
bulk Al to Al surfaces and surface steps, and to small 
Al molecules. The functionalities of the refined potential 
are found to extend significantly those of previous ones. 
As we are going to use an embedded atom interaction 
model, in this Section we briefly review the basic ideas 
of this approach, describe the details of the model, and 
finally assess its quality. 



A. Theory 

In the embedded atom method, each atom in a solid 
is viewed as an impxirity embedded in a host comprising 
all the other atomsO'EI The energy of the-host with im- 
purity is, according to Stott and Zaremba,£3 a functional 
of the unperturbed host electron density, and a function 
of the impurity type and position, 



E = Tz,r\ph(t)\, 



(1) 



where p/i(r) is the unperturbed host electron density, and 
Z and R are the type and position of the impurity. Here 
the energy of an impurity is determined by the electron 
density of the host before the impurity is added. The 
functional J- is a. universal fupction, independent of the 
host, but its form is unknown.til A simple approximation 
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to T is the so called local approximation, whereby the. 
impurity experiences a locally uniform electron density.c£l 
This can be viewed as the lowest-order term of an expan- 
sion involving the successive gradients of the density. The 
functional T is then approximated by 



E = F i { Pi (R i )) + ^^2<l>ij(Ri t 



(2) 



where fyj is a pair potential representing the electrostatic 
interaction, Rij is the distance between atoms i and j, 
and Fi denotes the embedding energy. The total energy 
of the system is a sum over all individual contributions: 



Eu 



= Y,F i (p h , i ) + ^^cl> ij (Ri j ). 



(3) 



A further simplification is introduced assuming that the 
host density ph^i at atom i is closely approximated by a 
sum of the atomic densities pj of the constituent atoms, 
i.e. p hti = J2j,(jjti) Pj(Rij)> witn Pj being the contribu- 
tion to the density at atom i from atom j. Equation 
(|J) is the form commonly used for molecular dynamics 
simulations of metals, and is known as embedded atom 
potential. 



B. Details of the Al-Al interaction potential 

The Ercolessi- Adams interaction model for Al was cop-, 
structed with the so called force matching methodtil 
and, in contrast to most other empirical models, it gives 
excellent structural and elastic properties for the bulk 
along with the correct surface interlayer relaxations at 
low-index surfaces. Furthermore, we found that the 
diffusion barriers for surface adatoms obtained by the 
Ercolessi- Adams_mpdel are in fair agreement with ab ini- 
tio calculations Ji3ilZl whereas those predicted by-, most 
other embedded atom potentials differ drasticallyES from 
ab initio results. We therefore started off from the 
Ercolessi- Adams potential to build our own refined Al- 
Al interaction. Without affecting the elastic properties 
and the surface relaxation properties, we introduced the 
following modifications to the model: 

f. An additional term was introduced in the pair 
potential ^ij m order to account for an expo- 
nential Bopi-Mayer-like repulsion at short Al-Al 
separation.113 This is a key requirement for studies 
of e.g. physical vapor deposition processes, where 
the energy of each single atom easily exceeds the 
thermal energy by as much as three orders of mag- 
nitude. 

2. In the low density region, three parameters of the 
embedding function F were changed in order to im- 
prove several reference quantities, namely the AI2 
binding energy, and vibrational frequency, and the 



adatom diffusion barrier height on the Al(lll) sur- 
face {i.e. the energy difference between the surface 
binding site and the saddle point). 

3. A fifth-order polynomial cut-off function was intro- 
duced, smoothly bringing the potential to zero at 
an interatomic distance of 5.56 A (slightly larger 
than the third-nearest neighbor distance in bulk 
Al). 

The total energy E tot of a system containing Al a toms 
in an arbitrary arrangement is written (see Sec. II A ) as 



E u 



The atomic density pi in arbitrary units is given as 
Pi = 2J P( r i]) x M r v,Ro,D ). 



(4) 



(5) 



The sum runs over all atoms that lie within the potential 
range i?o + -Do (5.56 A), which is enforced by the cutoff 
function / c (r, R, D). This function is zero for r exceeding 
R + D and unity for r less than R — D. For r within the 
interval (R — D,R + D) it is defined according to 



-5[4^ + l] 3 + l. 



(6) 



The function p(r) in Eq.(^) is spline- interpolated using 
the values reported in Table Q; the parameters Rq and 
Dq are given in Table The embedding function F(p) 
is also spline-interpolated, and the corresponding values 
for F(p) are collected in Table [II. 

The pair potential term in Eq.(ij) is written according 

to 

4>ij = [ <t>{nj) + (A exp{-Xnj} x f c {rij,It4,D^) 

-B)] x f c ( rij ,Ro,D ). (7) 



The function <f> is tabulated in Table IV. The first cut- 
off f c (rij , R(j>, D$) switches on the exponential repulsive 
term at small distances (r < 2.25 A), while / c (fy , i?o, -Do) 
terminates the interaction range of the potential. The 
corresponding parameters are given in Table ||. The ex- 
ponential term ensures that one gets a Born-Mayer repul- 
sion at short separations for, e.g., diatomic molecules.Ej 



C. Assessment of the potential 

We now present evidence that the potential just de- 
scribed yields satisfactory results for a variety of prop- 
erties of Al in different environments. In particular we 
address the bulk, the dimer, low-index surfaces and steps 
thereon, and self-diffusion on different low-index surfaces; 
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also included in this test section is the energy depen- 
dent sticking coefficient of high energy Al atoms on Al 
(111). The changes to the potential significantly im- 
proved agreement with experiment and other theoretical 
predictions in several instances where errors were typi- 
cally of order 50%. 



1. Bulk properties 



By construction our model does not alter the equilib- 
rium lattice constant ao, the cohesive energy E co h and the 
elastic properties of the previous model of Ercolessi. We 
obtain as in Ref. a = 4.03 A, E coh = 3.36 eV, C n = 
118 GPa, C12 = 62 GPa and C44 = 36 GPa. For compar- 
ison, the experimental valuesEll are Cn = 114 GPa i-CfSi 
= 62 GPa and C 44 = 32 GPa, and LDA calculationsO'Eil 
predict a = 3.98 A, E coh = 4.15 eV, Cn = 135 GPa, 
C12 = 70 GPa and C 44 = 35 GPa. 



2. Diffusion barriers 

Diffusion is central to many physical processes which 
determine the morphology of surfaces, such as step flow, 
nucleation, and growthHj It is of obvious importance to 
study diffusion processes theoretically, since direct ob- 
servations of surface diffusion by means of field ion mi- 
croscopy (FIM)E^I are limited to a few surfaces due to 
the respoase limits of the materials of interest to high 
voltages.c3 The barriers for single adatom difljasipH on 
Al surfaces calculated by Stumpf and Scheffleit3~tL3 us- 
ing ab initio LDA techniques, provide a stringent test 
for the present empirical Al-Al interaction model. It is 
generally accepted that diffusion on flat metal surfaces 
proceeds by either hopping or exchange. Ell In the two 
following subsections compare the results of the present 
model for these mechanisms with previous LDA results. 

Hopping Diffusion - During hopping diffusion the 
adatom is moving between minima of the potential en- 
ergy surface, i.e. between stable or metastable binding 
sites. On the (111) surface the stable adsorption sites 
are the 3-fold fee and hep sites; on the (100) surface 
there is single independent adsorption site, the four-fold 
hollow; the (110) surface is analogous, with a five- fold 
site. Hopping diffusion on the (110) surface is intrinsi- 
cally anisotropic, since it can proceed perpendicular or 
parallel to the [110] -oriented atomic rows, respectively 
via the short bridge or long bridge paths. The activa- 
tions energies for the long and short bridge are labelled 
E\\ and E±_ respectively. 

For each surface we performed total energy calcula- 
tions for the adatom sitting at the adsorption site and 
at the bridge site. At the latter site the total energy is 
minimized with respect to distance of the adatom from 
surface. All other Al positions are fully optimized. The 
energy difference between adsorption and bridge site is 



defined to be the activation energy for hopping diffusion. 
Further technicalities are discussed in the Appendix. Ta- 
ble summarizes the activation energies for surface self 
diffusion obtained with the present model and compares 
the barrier heights to results from previous ab initio cal- 
culations. The barrier for diffusion on Al (111) was used 
to set up the potential, hence it agrees with LDA results 
by construction. On Al (100), we obtain a diffusion bar- 
rier in good agreement with first principle calculations; 
we are not aware of experimental results for these self- 
diffusion barriers on Al(lll) and Al(100). According to a 
recent studyj23 the hopping-self-diffusion barriers calcu- 
lated ab initio in the generalized gradient approximation 
to density functional theory for unreconstructed fcc(100) 
surfaces equal one sixth the bulk cohesive energy; this 
is found to be the case for Al also in our calculations. 
In the case of Al (110), the present potential correctly 
predicts diffusion anisotropy, although some quantita- 
tive discrepanc»i-.esists with the-predictions of ab ini- 
tio calculations JiSlij FIM studiesa find no anisotropy for 
the diffusion on the (110) surface, and report a barrier 
height of 0.43 eV for both paths. This discrepancy with 
our result is due to the fact that here we only addressed 
hopping diffusion for demonstrative purposes. In fact, 
exchange diffusion normal to .the rows is known to have 
a barrier as low as ~ 0.6 eVjiaij which restores a rea- 
sonable agreement with experiment. 

Exchange Diffusion - Diffusion by atomic exchange oc- 
curs as the adatom replaces a surface atom, which in turn 
pops up at an adjacent stable surface site. Diffusien by 
exchange was discussed by Bassett and Webberc3 and 
Wrigley and EhrlichEH, and ,-aredicted theoretically for 
(100) surfaces by FeibelmanEa This diffusion mode can 
lower significantly the effective diffusion barriers. Views 
on why diffusion by exchange is favorable for some metal 
surfaces is discussed in Refs. |l7|, ^8|, ||c], and refer- 
ences therein. For Al (100), we calculated the activation 
energy for exchange diffusion (see Table ^) and found it 
lower than for hopping, as do first-principles calculations. 
The quantitative agreement is reasonably good. 



3. Surface and step energies 

The surface energy, defined as the difference between 
the energy of an atom at the surface and in the bulk 
environment, is usually calculated as 

E s \ah — N E'bulk 



77^arca 
^surf 



E. 



2A 

E s ub — N E'buik 



surf 



2A S , 



rf 



(8) 
(9) 



where -E s i a b is the total energy of the slab, N the total 
number of atoms in the slab, N sur { the number of atoms 
on each surface, -Ebuik is the total energy per bulk atom, 
A is the area of each free surface of the slab, and the 
factor 1/2 accounts for the two free surfaces of the sim- 
ulation cell; periodic boundary conditions are applied in 
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the planar directions. These formulas are not problem- 
free in general,Ej but we have checked that they are in 
the cases of interest to us. 

We calculated the formation energies for the low index 
(111), (100) and (110) Al surfaces with the present Al- 
Al interaction model. The comparison of our results to 
those obtained in ab initio LDA investigations, given in 
Table VI shows that the trends of surface energies of our 
Al model are consistent with the ab initio calculations. 
We obtain all surface energies about 20% lower than the 
LDA surface energies. Keeping in mind the known LDA 
overestimate of the blading energies (the LDA cohesive 
r 4 1 * 1 r T r L3 ' 3 about 20% larger than the ex- 

3V) the — * 
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{110}-facetted step 
on Al(100) 



{100}-facetted step 
on Al(111) 





{111}-facetted step 
on Al(100) 



{111}-facetted step 
on Al(111) 



FIG. 1. Atomic arrangement of the {110} and {111} 
facetted steps on Al (100), and the {100} and {111} facetted 
step on Al (111). 



corresponds to the height of the step. In the case of 
monoatomic-height steps, m is omitted in this notation. 

For the two low index surfaces Al (100) and Al (111), 
we calculated the formation energy of different steps. On 
Al (100) there exist two monoatomic steps, the close- 
packed {lll}-facetted and the more open {110}-facetted. 
The former belongs to the family of (1, 1, 2n + l)-surfaces, 
the latter to the family of the (1, 0, n)-surfaces.E3c£l 
For the calculation of the step energies we used the 
Al(l,l,15) = Al[9(100) x (111)] and the Al(l,0,9) = 
Al[9(100) x (110)] surfaces. 

On Al(lll) there are two types of close-packed steps, 
the {lll}-facetted and the {100}-facetted. The cor- 
responding vicinal surfaces belong to the (n, n, n — 2) 
and (n,n,n + 2) families, respectively. We used the 
Al(9,9,7) = Al[9(lll) x (111)] and the Al(8,8,10) = 
Al[9(lll) x (100)] surfaces. The geometry of the different 
steps on the AZ(100) and Al (111) surface is depicted in 
Fig.|l|. For all these vicinals, the terraces separating the 
steps have the same width of 9 atomic rows. We have ver- 
ified that step-step repulsion at these inter-step dists 
is already in the long-range elastic regime ~ d~^ s 
the steps are far enough to extract their formation energy 
without an unknown bias from the interstep interaction. 



Table VII lists the results for step formation energies, 
and compares them to ab initio data. The empirical Al 
potential describes the step formation energies for the 
two different steps on Al (111) in excellent agreement 
with first-principles calculations. More energy is needed 
to create steps on the close-packed Al (111) surface than 
on the more open Al (100). The open step on the Al 
(100) surface has a 20% larger formation energy than 
the close-p|SMsked step, in agreement with bond cutting 
arguments. 123 



4. The Al Dimer 



The energy per unit length of a step on a low- index sur- 
face is defined in terms of the energies of the low-index 
face and the vicinal surface used to simulate the step 
itself, and of the the geometrical parameters thereof:E3 

Estep — d s _ s -Evicinal ^terrace-^lc-w — index; (-^) 

with d s - s the step-step distance, and /terrace the terrace 
length on the vicinal surface. In analogy to the surface 
energy, the step energy can also be expressed per step 
atom. 

Although a stepped vicinal surface can be specified by 
its corresponding Miller indices, this notation is not very 
convenient, as it does not indicate at first sight the ge- 
ometrical structure of the surface. Thus we use instead, 
the notation [n(h, k, I) x m(h' , k', I')] by Lang et aL,EM-3 
where (h, k, I) and (h 1 , k', I') are the Miller indices of the 
terraces and ledges respectively; n gives the number of 
atomic rows in the terrace parallel to the step, and m 



The dimer is a stringent test for an Al-Al interaction 
model, since atoms in a dimer experience a very differ- 
ent chemical environment compared to bulk or surface 
atoms. For AI2, our model yields a binding energy per 
atom of 0.70 eV and a bond length of 2.70 A. The bind- 
ing energy was indeed used as input to determine the 
model parameters, matching the LDA valueO of 0.71 eV 
(-0.68 ± 0.03eV experimental, Ref. |oJ) The lowest vi- 
brational frequency is calculated to be v — 290 cm. -1 , 
in excellent agreement with the experimental valueEj of 
284.2 cm -1 . The predicted bond length of_.our model 
matches exactly the experimental estimate. EJ Thus, the 
present model describes satisfactorily the bonding of Al 
in the rather extreme case of the AI2 dimer. 

As the understanding of diffusion and growth requires 
a knowledge of the binding energies of small aggregates of 
adatoms, we also calculated the energy of two Al adatoms 
sitting at neighboring fee sites on an Al (111) surface. 
The energy gain with respect to isolated adatoms is 0.50 
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eV, i.e. the Al ad-dimer is stabilizedj-appreciably over 
separated adatoms. LDA calculationsEZI yield a similar 
energy gain of 0.58 eV. 



dent kinetic energy; the angle to the normal also has a 
drastic effect. Details of t he molecular dynamics calcula- 
tions are given in Sec. A 2. Further discussion and results 
on high energy deposition are reported in Rcf. [42. 



5. Sticking coefficient for hyperthermal Al atoms 

During physical vapor deposition the Al atoms emitted 
from the sputter source have a non-thermal energy dis- 
tribution, with kinetic energies exceeding 10 eV. There- 
fore the sticking coefficient, a key ingredient for a reliable 
modelling of metal film growth, cannot be assumed to be 
constant and independent of the particle's energy as it is 
typically done. In order to elucidate the dependence of 
the sticking coefficient on impingement energy, we start 
our simulations with the incident Al atom placed outside 
the interaction range of the surface. Its initial kinetic 
energy is set in the range of to 125 eV, and its starting 
angle off the surface normal in the range 0° to 60° , which 
corresponds to typical ionized physical vapor deposition 
conditions. The trajectories of the incident atom, and 
of any other atom which may be etched away from the 
surface upon impact, are monitored until either a certain 
time span has elapsed, or the outcoming atoms (in the 
case of reflection or etching) have traveled a distance of 
10 A away from the surface. Analyzing 200 trajectories 
per incident energy and angle, we collected a statistically 
significant sample of well-defined adsorption, reflection, 
and etching events. The relative probability of the stick- 
ing coefficient is calculated as the ratio of the number of 
adsorption events to the total number. The typical sta- 
tistical error in the reaction probability thus determined 
is below 5 %. Fig.|^ depicts the sticking coefficient as a 
function of energy for Al atoms imping normally on the 
surface (solid circles) or at an off-normal angle of 40° 
(open circles). 
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FIG. 2. Sticking probabilities for hyperthermal Al atoms 
impinging on an Al (111) surface. As a function of kinetic en- 
ergy the solid line with filled circles depicts the reaction prob- 
abilities for Al atoms impinging normal on the surface and the 
dashed line with the open circles show the latter quantity for 
an angle of 40° to the normal. 

The sticking probability varies strongly with the inci- 



III. APPLICATION: FREE ENERGY 
CALCULATIONS 

In this applicative Section of the paper, we compare 
surface free energies of Al computed using different levels 
of quasi-harmonic approximation, and thermodynamic 
integration via Monte Carlo simulations. The latter ef- 
fectively functions as "exact" reference for the various 
harmonic approximation. Before presenting the results, 
we briefly review the background theory of the different 
approaches. 

A. Theory 

Thermodynamic Integration - The free energy can- 
not be calculated as an ensembie,|avefage.c3 The method 
of thermodynamic integration! 1 ^rTcZI circumvents this 
problem starting from the concept of Stockmayer fluid, 
a fictitious system in which the interparticle interaction 
potential U\ is gradually switched on from a known refer- 
ence potential U% to the actual, full interaction potential 
U; the mixing of Ut and U into the effective potential is 
controlled by a parameter A: 

U x = (l- X)U h + XU. (11) 

The key relation of the method concerns the derivative 
of F with respect to A: 

The subscript A means that the average has to be eval- 
uated with the interaction potential U\. Integrating the 
latter equation one arrives at the following expression for 
the free energy of the system of interest: 

F x =i = F x=0 + f (U- U h ) x d\ (13) 
Jo 

The usual choice for the reference system is the Einstein 
crystal (i.e. a system of non-interacting harmonic oscilla- 
tors with the interaction potential Uh — \muir) |r< — 
r i0 | 2 ), whose free energy is 

F A=0 = -3iYfc B Tln(J-), (14) 

with @ D the Debye temperature (394 K for A@). The 
free energy of the real system can thus be obtained at any 
given temperature by a series of canonical Monte Carlo 
simulations. 
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A faster way to obtain the temperature variation of the 
free energy is to integrate the thermodynamical relation 



d_ (F_ 



\rp2l1 



(15) 



from a reference temperature upwards. This requires a 
simple {e.g.) Monte Carlo ensemble average of the energy 
for each temperature, and of course a reference value of 
F from thermodynamic integration. 

In summary, the free energy of the system at a refer- 
ence temperature To is determined with Eq. (fl3|) ; then, 
the temperature variation from To to T is calculated us- 
ing Eq. (|l5|) . Both steps were performed by canonical 
Metropolis Monte Carlo simulations.e3E3 As detailed be- 
low, thermal expansion is taken into account performing 
the NVT Monte Carlo calculations at the temperature- 
dependent lattice constant determined independently by 
NPT molecular dynamics. 

Quasi-harmonic approach - A popular approach to free 
energy calculations is the quasi-harmonic approximation. 
Thereby, the full interatomic potential is replaced by its 
quadratic expansion about the atomic equilibrium po- 
sitions. The system is then equivalent to a collection 
of harmonic oscillators, and diagonalization of the cor- 
responding dynamical matrix yields the squares of the 
normal- mode frequencies, i.e. the phonon spectrum. In 
bulk systems the dynamical matrix is a 3 x 3 matrix; for 
a slab system it is a 3£ x 3£ matrix, where £ is the num- 
ber of atomic layers in the slab. The dynamical matrix 
is given bytj 



1 E 

m. z — ' 



$ Q/3 («')exp[ iq (r -r )] (16) 



where the force constant matrix $ a p(£ £') is defined as 



d 2 u 



du a {£) dup{£') 



(17) 



The subscript "0" indicates that the second derivatives 
are to evaluated at the true mean positions of the atoms, 
with any displacements from the bulk positions {e.g. sur- 
face relaxations or reconstructions) taken into account. 
The equilibrium positions Xq,i/q, Zq of the atoms are given 
by the vectors I — {l\,li,lz)\ the £3 axis is perpendic- 
ular to the surface and the position of an atom within 
a plane is specified by Si, £2, and u a {£) describes the a 
component {a = x, y, z) of the position of the ^-th atom 
from its mean position Xq, j/q, Zq. 

The phonon spectrum of bulk Al and of Al (100) are 
displayed in Fig.|3j, upper and lower panels respectively 
(for both calculations supercells comprising 50 layers 
stacked along (100) have been employed). A variety of 
surface modes appear in bulk gaps or split off from bulk 
band edges. These additional modes are the source of 
the different vibrational free energy of surface systems 
in comparison to bulk systems. The free energy in this 



approximation is calculated for lattice and geometrical 
parameters a at temperature T as 



F(a, T) = Eo{a) 



fc B T^ln(2 sinh 



2k R T 



). (18) 



The sum runs over all phonon polarizations j and wave 
vectors k in the Brillouin zone, with u>j(k) the frequency 
of the corresponding modes. Both the frequencies, and 
the internal energy Eq(sl) of the ideal static lattice, de- 
pend on all the lattice and geometrical parameters a. 
The latter include the bulk lattice constant and, for the 
surface, the additional geometrical parameters involved 
in relaxations or reconstructions. 
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FIG. 3. Bulk (a) and surface (b) phonon spectra for Al 
calculated with a 50 layer slab. In (b) the slab has two (100) 
surfaces, the frequencies are plotted along lines of high sym- 
metry. The corresponding 2d-Brillouin zone is shown in the 
inset of panel (a). 

In a bulk system the forces on each atom are zero by 
symmetry, independently of a, so that it is strictly correct 
to neglect the first derivatives in the quadratic expansion 
of the potential energy. We have verified that the quasi- 
harmonic approximation does indeed work very well for 
the bulk even in comparison to thermodynamic integra- 
tion. For a surface, the situation is different, since the 
interlayer spacings (expecially those of the top surface 
layers) will change with respect to the bulk, i.e. the sur- 
face will generally either contract or expand. The aver- 
age equilibrium positions of the near-surface layers (the 
interlayer spacings in the case of simple relaxation) are 
determined by the minimum of the free energy, but that 
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is, of course, unknown a priori. In addition, the harmonic 
expansion is not strictly correct (since the forces, i.e. the 
derivatives of the internal potential energy are not zero 
at the free energy minimum), which is why one calls this 
the quasi-harmonic approximation to begin with. In fact, 
it is clear that there are several levels of approximation 
possible for the quasi-harmonic approach; here we con- 
sider some of those: 

1. The computationally simplest way is to optimize 
the atomic configuration and calculate the phonons 
at zero temperature, and evaluate the free energy 
vs. T using those ingredients for all T. Within 
this approach, the quasi-harmonic approximation 
is strictly valid, as we expand the potential energy 
function around the equilibrium positions, and the 
T dependence enters solely with Eq.([l8|). In real 
systems, of course, both the surface internal en- 
ergy and the vibrational contribution to the free 
energy will change with T, but it is a priori unclear 
to what degree this influences the result. 

2. Another way to account for the effects of finite tem- 
perature is to take the T = atomic positions, 
rescale their coordinates as a function of the tem- 
perature according to the appropriate bulk thermal 
expansion coefficient, and recalculate the phonons 
(and hence the free energy) for the expanded lat- 
tice. This is a hybrid case in which T not only 
affects the force constants, but also the surface in- 
ternal energy. Of course it is arbitrary to use the 
scaled T=0 interlayer spacings at non-zero tem- 
peratures. Also, it should be kept in mind that 
the harmonic approximation is not strictly valid 
for the expansion of the potential energy around 
non-equilibrium positions. 

3. A further possibility is to rescale all the coordi- 
nates according to thermal expansion first, and 
then re-optimize all atomic positions; the phonons 
and the free energy are calculated for that geom- 
etry. Here one is consistent with the prerequisites 
of the quasi-harmonic expansion, but at the cost of 
getting wrong interlayer spacings at the surface. 

4. The real thing is of course to minimize the total 
free energy with a "self-consistent" adjustment of 
the atomic positions of all layers in the slab system, 
resulting in the thermodynamic equilibrium con- 
figuration of the surface system. In practice, one 
starts with the bulk positions rescaled according 
to thermal expansion, and then adjusts the inter- 
layer spacings of a few near-surface layers to obtain 
the minimum of the free energy. The major con- 
tribution is generally due to the first two surface 
layers, the only having a sizable displacement from 
the bulk interlayer spacing. In our calculations, we 
therefore changed d\2 by ±3% and d^3, by ±2%. 



B. Results 



We now compare the free energy of an Al(100) sur- 
face calculated within the different quasi-harmonic ap- 
proaches (1-4) described above, with the results of ther- 
modynamic integration; the latter effectively functions as 
exact reference since it takes the full potential into ac- 
count, hence in particular all anharmonic contributions. 
We chose the (100) surface for demonstrative purposes, 
as intermediate between the closed packed (111) and the 
more open (110) surface. 



1. Comparison of different methods 



In order to obtain the bulk lattice constant at different 
temperatures we first nfiaJhrmed zero-pressure molecular 
dynamics simulatiorO'Ej'Ej using our Al potential. In the 
linear regime the expansion coefficient is a = 1.64 x 10z£ 
A/K, the experimental value being 2.36 x 10~ 5 A/K.Ej 
Deviations from linearity set irO at about 500 K. The 
dimensions of the simulation cell with periodic bound- 
ary conditions correspond, in all subsequent calculations, 
to the bulk lattice constant at the relevant temperature. 
The surface free energy calculations imply the evaluation 
of the bulk free energy, and of the free energy of a slab 
system with two surfaces. The surface free energy is then 
determined with Eq.(^|). 

In FigJJ we compare the surface free energy calcu- 
lated with thermodynamic integration and the quasi- 
harmonic approach. All versions of the latter underes- 
timate severely the temperature variation of the surface 
free energy. This is mainly due to the neglect of anhar- 
monicity, which is also responsible for thermal expansion. 
A first important result is then that at temperatures 
T > 0_d, the harmonic approximation is inappropriate 
for Al surfaces. 




500 600 700 800 900 
T[K] 
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FIG. 4. Surface free energies for Al (100) calculated with 
different quasi-harmonic approaches (as discussed in the text) 
and with the method of thermodynamic integration. Solid 
line: approach (1), zero-temperature phonons for all temper- 
atures; dotted line: approach (2), positions rescaled according 
to thermal expansion; dash-dotted line: approach (3), as (2) 
with re-optimization of atomic positions; crosses: approach 
(4), minimization of the free energy in the {di2,d23} plane; 
dashed line: (TI) thermodynamic integration, reference for 
harmonic approximations. See text for more details. 



Note that the failure of the harmonic approximation 
for the present relatively high-temperature calculations 
does not affect the successes of this approach at low tem- 
peratures, an example being, the recent first-principles 
calculations for Be surfaces£3 The reason why those re- 
sults are compatible with ours is clearly that we work 
well above the Debye temperature of our system (~ 400 
K), whereas the highest temperature considered in Ref. 
|| is 750 K, well below ~ 1000 K (as extracted from 
a Dcbye-Einstein model). Of course, the quasi-harmonic 
approach will generally fail if applied to systems at suffi- 
ciently high temperatures. 




FIG. 5. Panel (a): dependence of the plain surface energy 
on the interlayer spacing di2 and ^23. Panel (b): vibrational 
contribution to the surface free energy as a function of the 
interlayer spacings. Panel (c): total surface free energy, sum 
of the two previous contributions. 



To sort out the relative merits of the various levels of 
harmonic approximation, we focus on the effects of the 
interlayer spacing d\2 (between first and second layer) 
and c?23 (between second and third layer) on the surface 
energy and on the vibrational contribution to the sur- 
face free energy. Essentially this is the fourth level of 
approximation mentioned earlier. We pick T=450 K for 
demonstrative purposes, and expand the lattice accord- 
ingly. 



Panel (a) in Fig.|| shows the variation of the plain 
surface energy as a function of the interlayer spacings 
(expressed in turn in percentage of the bulk interlayer 
spacing). An increase of the interlayer spacings tends to 
increase the energy drastically, a decrease to reduce it. 
The minimum at around -3% for both spacings. These 
fairly unrealistic values result from the optimization of 
the interlayer spacing for a laterally expanded surface. 



The excess surface free energy, i. e. the vibrational con- 
tribution is shown in panel (b) of Fig.|^. An increase 
of the interlayer spacings leads to softer force constants 
and hence to lower frequencies, which yield according to 
Eq.(|l8|) a more negative value of the surface excess free 
energy. The dependence on the first spacing is stronger 
than on the second, although conceptually both spacings 
would tend to positive infinity (decoupled Al planes) if 
only the vibrational contribution mattered. 

The opposing tendencies of the internal and vibrational 
contributions tend to compensate; in fact, summing up 
the plain surface energy and the vibrational contribution, 
one arrives at the total surface free energy depicted as a 
function of the interlayer spacings in panel (c) of Figj^. 
For the Al (100) surface, the minimum in the free energy 
corresponds to d 12 = —0.5 % and ^23 = —0.9%, a com- 
promise between the gain of free energy upon outward 
relaxation, and that in plain surface energy upon inward 
relaxation. 

In approaches (1) and (2) from our above list 
(force constants from zero-temperature or rescaled zero- 
temperature posizions), the interlayer spacings are d\2=- 
1.5% and (i23=-1.3 %. These are rather close to the min- 
imum of the free energy found by direct minimization in 
approach (4); indeed, with reference to Fig.[|, both ap- 
proach (1) (solid line) and (2) (dashed line) match rather 
closely the values of approach (4) (crosses) , i. e. of the full 
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quasi-harmonic calculation. Approach (3) (dash-dotted 
line), where we rescaled the lattice constant and then re- 
optimized all atomic positions, fails badly, going astray 
already near the Debye temperature, and progressively 
more so for higher temperatures. This is due to the in- 
correct (free-energy-wise) spacings imposed on the near- 
surface layers by the minimization of the internal energy. 
The spacings are found to be d\2=— 3.4% and d23=~3.1 
% at 450 K, and d 12 =-5.8% and d 23 =-5.5 % at 900 K. 
A glance at panel (c) of Fig.|] reveals that both of these 
points in the {di2, ^23} plane do indeed correspond to free 
energies very far away from the minimum (especially at 
the higher temperature). 

In conclusion, the most naive and simplest approach 
of exporting the T = force constants and surface en- 
ergy to non-zero temperature does indeed underestimate 
considerably the temperature variation of the surface free 
energy with respect to thermodynamic integration, but 
it gives gives an agreement comparable to, or better than 
the sophisticated adjustment of the interlayers to find the 
free energy minimum. 



slabs with periodic boundary conditions for the lateral 
cells. The supercells contained of 672, 550 and 560 atoms 
for the Al(lll), Al(100) and Al(110) surfaces and con- 
sisted of 12, 11 and 9 atomic layers. In order to determine 
the surface energies of Al(lll), Al(100) and Al(110) the 
supercells contained 1080, 550 and 560 atoms, arranged 
in 9, 11 and 16 atomic layers. The step formation ener- 
gies were obtained from systems containing 4 steps and 
72, 105, 102 and 102 atoms per layer corresponding to a 
total number of 1312, 2724, 1368 and 2532 atoms for the 
Al(l, 0, 9), Al(l, 1, 15), Al(8, 8, 10) and Al(9, 9, 7) surface. 
All forces F per atom have been brought below a thresh- 
old of 10 -5 eV/A. We estimated the errors in the total 
energies due to the finite supercell size to be well below 
10" 4 eV/atom. 



2. Molecular dynamics calculation of the sticking 
probability 



IV. SUMMARY 

We have presented a refined Al interatomic poten- 
tial for classical dynamics and Monte Carlo simulations. 
We thoroughly tested its functionalities, finding it to be 
very accurate for a variety of systems. Next, we applied 
it to evaluating the performance of quasi-harmonic ap- 
proaches to free energy calculations for surfaces, compar- 
ing the latter results with full thermodynamic integration 
results. For Al surfaces, the quasi-harmonic approxima- 
tion shows a progressively increasing error for tempera- 
tures above <3d- Different levels of quasi- harmonic ap- 
proximation have been compared; for Al, the simplest 
method of using zero-temperature phonons to compute 
the free energy at all temperatures is as accurate as the 
explicit minimization of the free energy with respect to 
geometrical parameters. 
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APPENDIX A: COMPUTATIONAL DETAILS 

1. Substrate sizes 

For the calculation of the surface self-diffusion barriers, 
and surface and step energies, we have employed finite 



The reaction probabilities were calculated in classical 
molecular dynamics simulations using our Al interaction 
potential. The integration was performed with a 5-th or- 
der Runge Kutta method with an adaptive timestep, in 
order to ensure total energy conservation throughout the 
simulation. Supercells containing 1320 atoms arranged 
in 10 atomic layers were employed; cell dimensions are 
chosen so as to avoid artifacts of the in-plane periodicity. 
The starting configuration is chosen to be a (111) surface, 
the one Al surface with the lowest formation energy. All 
atomic coordinates are allowed to evolve dynamically, ex- 
cept those of the two bottom layers of the supercell. The 
surface temperature is set at 450 K (i.e about 1/2 of the 
melting temperature, and ~ 15 % larger than the bulk 
Debye temperature). 



3. Monte Carlo calculations within the canonical 
ensemble 



All Monte Carlo calculations were be performed within 
the canorueaL,ensemble, using the standard Metropolis 
technique J13'c2l The maximum atomic displacement was 
automatically adjusted in order to get an acceptance ra- 
tio of 0.4. It was not systematically studied that this 
acceptance ratio was an optimum, but well converged 
statistical averages were obtained with a typical number 
of Monte Carlo moves of order 10 4 times the number of 
atoms in the system. Before averaging, the system was 
equilibrated for a number of steps of order 500 times the 
number of atoms in the system. For the Al (100) surface 
we used in total 384 atoms. 







4. Quasi-Harmonic free energy calculations 



Within the quasi-harmonic methods we employed slab 
geometries with 20 atomic layers each containing 32 
atoms. For the k-space summation we used grids typ- 
ically containing 2500 equally spaced k-points. Careful 
tests showed that this number of k-points yields well con- 
verged results. 
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TABLE I. Parameters used to define the atomic density 
function p(r). The positions of the spline knots and the val- 
ues at the knots are given. Also the first derivatives at the 
first and last knots are given. 
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r [A] PW P'(r) [1/A] 

0.0000 0.000 0.0000 

1.8000 6.3820X10" 1 
1.9000 7.6541 xlO" 1 
2.0211 8.6567xl0 _1 
2.2737 9.2521 xlO -1 
2.5264 8.6200xl0 _1 
2.7790 7.6273X10- 1 
3.0317 6.0648X10" 1 
3.2843 4.6603xl0 -1 
3.5370 3.3874X10" 1 
3.7896 2.3257xl0 _1 
4.0422 1.0905 xlO -1 
4.2949 5.2491 xlO" 2 
4.5475 3.9170xl0 -2 

4.8001 3.0828xl0~ 2 
5.0528 2.5021 xlO -2 
5.3054 1.4722 xlO" 2 

5.5600 0.0000 2.1298x10" 



TABLE IV. Parameters used to define the pair potential 
4>{r). The positions of the spline knots and the values at the 
knots are given. Also the first derivatives at the first and last 
knots are given. 



r [A] cj>(r) [eV] g/(r) [eV/A] 

2.0211 1.9601 -7.0273 

2.2737 6.8272X10- 1 

2.5263 1.4737xl0 _1 

2.7790 -1.8818xl0~ 2 

3.0317 -5.7601xl0~ 2 

3.2843 -5.1984xl0" 2 

3.5369 -3.7635xl0 -2 

3.7896 -3.7373xl0 -2 

4.0422 -5.3135xl0~ 2 

4.2949 -6.3286xl0~ 2 

4.5475 -5.4810xl0~ 2 

4.8001 -3.7288xl0" 2 

5.0528 -1.8887xl0 -2 

5.3054 -5.8523xl0 -3 

5.5600 0.0000 5.9065xl0 -6 



TABLE II. Parameters entering the potential energy func- 
tion from Eq.(Q) 

Parameter 



Ro [A] 
Do [A] 
R* [A] 
D* [A] 
A [eV] 
A [1/A] 
B [eV] 



5.46 
0.10 
2.00 
0.25 
7255.44 
4.42085 
1.04897 



TABLE III. Parameters used to define the embedding 
function F(p). The positions of the spline knots and the val- 
ues at the knots are given. Also the first derivatives at the 
first and last knots are given. 



P 



F(p) [eV] 



u'(p) [1/eV] 



0.0 
0.1 
0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
1.0 
1.1 
1.2 
1.4 



0.0000 
-0.8139 
-1.2697 
-1.6799 
-2.0296 
-2.2520 
-2.4272 
-2.5517 
-2.6052 
-2.6440 
-2.6571 
-2.6456 
-2.6087 
-2.4525 



-12.375 



1.0620 



TABLE V. Comparison of selected hopping and exchange 
diffusion barriers on low-index Al surfaces obtained with the 
present model and in ab initio calculations. Al (111) is also 
included for completeness. 



This work 



ab initio 



Al (111) hopping 


0.04 


0.04 a 


Al (100) hopping 


0.60 


0.68 a ,0.65 b 


Al (100) exchange 


0.50 


0.35 a 


Al (110) 1 hopping 


1.13 


1.06 a 


Al (110) || hopping 


0.30 


0.60 a 


a Reference b Reference ^8] 




TABLE VI. Surface 


formation energies for low index Al 


surfaces calculated with the present model and with other 


theories. 






This work 


LDA a 


System (eV/atom) 


(eV/A 2 ) 


(eV/atom) (eV/A 2 ) 


Al(lll) 0.38 


0.054 


0.48 0.070 


Al(100) 0.48 


0.059 


0.56 0.071 


Al(110) 0.74 


0.065 


0.89 0.080 



1 Reference |l7| 



TABLE VII. Step formation energies for low index Al sur- 
faces calculated with the present model and with other theo- 
ries. 



This work 
System (eV/A) (eV/atom) 



LDA a 
(eV/A) (eV/atom) 



Al[9(100) x (111)] 0.055 0.142 

Al[9(100) x (110)] 0.066 0.240 

Al[9(lll) x (111)] 0.083 0.215 

Al[9(lll) x (100)] 0.085 0.222 j 



0.082 
0.088 



0.232 
0.248 



1 Reference |l6| 
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temperature 



expected 



calculated 




*■ temperature 




